get.change.index <- function(x){	

  ## Number of observations
  m <- length(x)

  ## Mean value of the vector
  xbar <- mean(x)
  ## Summing functioning
  s <- rep(0,m+1)

  for(i in 1:m){
    s[i+1] <- s[i] + x[i] - xbar
  }

  smax <- max(s) - min(s)
  s <- abs(s)

  index <- which(s == max(s))
  r <- data.frame(change.index = index, change.statistic = smax)

  return(r)

}

get.change.confidence <- function(x, return.data = FALSE, return.is.change = FALSE, ...){
	
	## Check if the confidence level was specified
	change.confidence <- match.call(expand.dots=TRUE)$change.confidence
 		
	if(!is.numeric(change.confidence)){
		change.confidence <- 0.95
	} else if(change.confidence >= 1) {
		change.confidence <- 0.99
	} else if(change.confidence < 0.50) {
		change.confidence <- 0.50
	}
	
	
	## Check whether the number of iterations size was specifiied
	n <- match.call(expand.dots=TRUE)$n

	if(!is.numeric(n)){
		n <- 1000
	} else if(n < 10) {
		n <- 1000
	}
	

	## Create the boostrap and determine whether change occurred
	## the r vector stores the change test statistic for the n samples
	r <- rep(0,n)
	for(i in 1:n){
		r[i] <- get.change.index(sample(x, replace=FALSE))$change.statistic
	}
	
	## quantify the change confidence
	## the test statistic is the number of times the resample data has an smax > the observed data smax
	
	observed.data.change.index <- get.change.index(x) 
	change.confidence.observed <- 1-length(r[r > observed.data.change.index$change.statistic])/length(r)
	
	## 3 options to return the 
	if(return.data){
		## r holds the test statistic for each resample
		return(r)
	} else if(return.is.change) {
		if(change.confidence.observed >=change.confidence){
		  df <- data.frame(observed.confidence = change.confidence.observed, is.change = TRUE) 
		} else {
		  df <- data.frame(observed.confidence = change.confidence.observed, is.change = FALSE) 
		}
		return(df)
	} else {
		return(change.confidence.observed)
	}
	
}

plot.change.confidence <- function(x){

  t.baseline <- get.change.index(x)$change.statistic
  t.bootstrap <- get.change.confidence(x, return.data=TRUE)
  mtxt <- "Histogram of CUSUM Test Statistic"
  h <- hist(t.bootstrap, xlim=c(min(t.baseline,t.bootstrap),max(t.baseline,t.bootstrap)), main=mtxt, col="blue")
  points(c(t.baseline,t.baseline),c(0, max(h$counts)), type='l', lwd=2, col='red')

}

get.change.interval <- function(x, change.index){
	
	## x: observations
	## change.index: location within x where change occurred
	indexA <- 1:(change.index-1)
	indexB <- (change.index):(length(x))

	n <- 1000
	ci.count <- 0
	mean1 <- mean(x[indexA])
	mean2 <- mean(x[indexB])
	error <-  c(mean1 - x[indexA], mean2 - x[indexB])
	
	change.point.location <- rep(NA, n)
	for(i in 1:n){
		x1 <- mean1 + sample(error, length(indexA), replace = TRUE)
		x2 <- mean2 + sample(error, length(indexB), replace = TRUE)
		x1x2 <- c(x1, x2)
		change.point.location[i] <- get.change.index(x1x2)$change.index
	}
	
	return(change.point.location)

}
